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We present in this paper the numerical treatment of the coupling between hydrody- 
namics and radiative transfer. The fluid is modeled by classical conservation laws (mass, 
momentum and energy) and the radiation by the grey moment Mi system || . The scheme 
introduced is able to compute accurate numerical solution over a broad class of regimes 
from the transport to the diffusive limits. We propose an asymptotic preserving modifica- 
tion of the HLLE scheme in order to treat correctly the diffusion limit. Several numerical 
results are presented, which show that this approach is robust and have the correct behav- 
ior in both the diffusive and free-streaming limits. In the last numerical example we test 
this approach on a complex physical case by considering the collapse of a gas cloud lead- 
ing to a proto-stellar structure which, among other features, exhibits very steep opacity 
gradients. 

Key Words: Radiative transfer; radiation- hydrodynamics; diffusion limit; 
transport; hyperbolic system; Godunov-type methods 

1. INTRODUCTION 

Radiation hydrodynamics plays an important role in astrophysics, laser fusion 
and plasma physics. For many years, efforts have been underway to develop mathe- 
matical models and numerical schemes to obtain accurate predictions at reasonable 
computing cost in this domain. 

One of the main difficulty is to obtain accurate numerical computations in the 
various regimes that can be encountered due to values of material opacities which 
can vary from several order of magnitude. This can be achieved by using the full 
radiative transfer equation but it is still out of range for complex simulations in 
particular in two or three spatial dimensions. To overcome this difficulty several 
models have been derived. For large values of the material opacity, an asymp- 
totic analysis leads to hyperbolic/parabolic systems of equation refered to as the 
equilibrium-diffusion limit or non equilibrium diffusion approximation ( p2| ,]jl|, 
JlO[ , 0]). On the other hand, for smaller values of the opacity the streaming limit 
can be handled by using an hyperbolic system of equations coupling some closure 
of the moment system of the transfer equation with the hydrodynamic system (Jl2|, 
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PI , pH,[p]) . Unfortunately for several real applications, such as radiative shocks 
or star formation, very different regimes are encountered in the same simulation. 
Coupling different models on various zones introduces large drawbacks due to the 
domain partition and some loss of accuracy in the transitional zone. 

From a numerical point of view high-resolution schemes are necessary to take 
into account the different space and time scales. Godunov-type schemes based on 
exact or approximate Riemann solvers are efficient for shock problems (Q, 0], 
[fl3|| ). In radiation hydrodynamics, without source terms, hydrodynamic and radi- 
ation equations decouple so that a Riemann solver for the entire system can easily 
be derived from Riemann solvers for the hydrodynamic subsystem and for the ra- 
diation subsystem. Nevertheless, the difficulty of this approach is to obtain the 
correct diffusion limit for the scheme and correct behavior in transitional regime. 
Let us notice that the regime at the numerical level depends on the ratio of the 
mean free path of photons and the size of the local cell. An other difficulty is 
related to the time scales in the considered problems. In most cases radiation has 
to be treated implicitly since the time step given by usual CFL condition is much 
too small. When possible the hydrodynamics may be treated explicitly (||), but 
in some problems the time step given by the CFL condition on the hydrodynamic 
subsystem is still too restrictive and hence a fully implicit approach is necessary. 

In this paper we propose new model and numerical scheme to describe radiation 
hydrodynamics in a wide range of regimes. The equations for the fluid are the clas- 
sical conservation laws (mass, momentum and energy). The radiation is described 
by the grey moment system of the transfer equation (radiative energy and radiative 
flux) with the Mi closure introduced in (|5|). This closure, based on the minimum 
entropy principle, is able to describe large anisotropy of the radiation while giving 
also the correct diffusion limit. The Eddington factor is a non constant function 
of the reduce flux and therefore the radiation subsystem is a nonlinear hyperbolic 
system. Here, these equations are written in the comoving frame, introducing some 
non conservative coupling terms. A Riemann solver is developped by performing a 
wave decomposition neglecting the nonconservative coupling terms and the source 
terms, which is trivial since the two subsystems then decouple. Afterwards the 
nonconservative terms are reintroduced in the scheme, by using in the radiation 
equations the velocity given by the hydrodynamic solver. This approach, natural 
in the streaming regime, can be improved to capture the diffusion limit. For that 
purpose the Riemann solver for the M\ subsystem is modified to take into account 
the stiff relaxation term along the ideas developped in (||). With this modifica- 
tion our approach is able to correctly treat the streaming regime as well as the 
diffusion regime with the same solver. Finally as radiative shocks or star formation 
simulations need high resolution and very different time scales we describe the im- 
plementation of our radiation hydrodynamics solver in a moving grid-time implicit 
framework. 

The outline of this paper is as follows. In the second section the equations of 
radiation hydrodynamics with the Mi-closure are given. In section 3 we present 
the numerical methods for solving these equations. The actual implementation is 
detailed in section 4. Numerical experiments are shown in section 5, first on radia- 
tive shocks to demonstrate the validity of our method in various regime, then on 
the collapse of a gas cloud leading to a stellar-like structure. The last section is the 
conclusion of the paper. 
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2. THE PHYSICAL MODEL 



The equations of grey radiative hydrodynamics under the assumption of local 
thermodynamic equilibrium can be found in ([^H). We have choosen in this work 
to use the radiative equation in the fluid (or comoving) frame : 



d tP + V-H = (1) 

d t pu + V.Lou® u + PI] = — pF r + F (2) 

c 

d t E + V.[u(E + P)} =-K P c{a r T i -E r ) + {-pF r + F).u (3) 

c 

d t E r + V.[E r u] + W.F r + P r : Vu = K P c(a r T 4 - E r ) (4) 

d t F r + V.[F r u] + c 2 V.P r + (F r .V)u = -npcF r (5) 

where p is the matter density, u the velocity, E the total matter energy density, P 
and T the pressure and the temperature of the material, F the force density, E r 
the radiative energy density, F r the radiative flux, P r the radiative pressure tensor 
and k the mean grey opacity. 

The first part of the system ([|)-(|]) describe the evolution of the matter and the 
second one (f|)-(||) the evolution of the radiation. These two fluids are coupled by 
the exchange of energy and momentum through respectively the Kpc(a r T 4 — E. r ) 
and K P F r /c terms. 

In order to close system ([!])- (||) it is necessary to give two equations of state for 
the flow and for the radiation. For the flow fields we assume a 7-law p — (7 — l)pe 
(with 7 = 5/3 in our applications). For the radiation fields we need a closure 
relation giving the radiative pressure P r in term of the radiative energy E r and 
of the radiative flux F r . A rather simple closure is to assume that the radiative 
flux is isotrope which leads to the closure P r = \E r . Such an approximation is by 
construction not very well suited to model flow containing large anisotropy in the 
radiative flux which is generally the case for radiative shocks. We will therefore 
use the closure given by the Ml model j5j. In this model, the radiative pressure is 
given by: 



P r = T>E r (6) 
where the Eddington tensor D is defined by: 

T>= l —>Ll + i -*-±n®n (7) 

X = ^ffl (8) 
5 + 2V4-3I/I 2 

I is the identity matrix, \ t ne Eddington factor, / = the reduce flux and 
n = //|/| is a unit vector aligned with the radiative flux. 
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Let us remark that we can write the conservation of total (i.e. radiative and 
material) momentum and energy. Neglecting terms in (u/c) 2 and terms involving 
the fluid acceleration this leads to the following equations 



d t (pu+^F r ) + V.[pu®u + PI + P r +u^-} = F-(J.V)u (9) 

dt(E + E r +u^)+W. [u{E + E r + P + P r ) + F r ] = F.u (10) 

Some authors make the choice to discretize directly the total mass, momentum 
and energy conservation equations with either the material or radiative momentum 
and energy equations to close the system. In this paper we prefere to work with the 
system (|l|)-(^) where the flow and the radiation play a more symmetric role. The 
coupling between the flow and the radiation is done by the relaxation source terms 
and by the terms introduced by the comoving frame. The mathematical analysis 
of the hyperbolicity of this system do not take into account the source terms. The 
Jacobian A of the system (0)-@ can be written 

_4 = (Ahh Arj 
y-A r / ? A rr J 

where the blocks are associated with the partitioning of unknowns between hydrody- 
namic variables (denoted h) and radiative variables (denoted r) . The hydrodynamic 
diagonal block Ahh is classical and writes 

/ 1 \ 

A hh = (7-3)£ -(7-3)« 7-1 (12) 

\-u(h-( 7 -1)^ H-(j-1)u 2 U1 J 

where H = e + p/p + u 2 /2 is the total enthalpy of the flow. Since the radiative 
variables do not appear in the left-hand-side of the hydrodynamic equations we 
have 

/0 0\ 

A hr = , (13) 
\o Oj 

The radiative diagonal block A rr for the Ml subsystem in the comoving frame 
depends on the velocity of the flow and writes 

Ar= (c 2 (x(/)"X'(/)/) u + cx'U))- (U) 

Finally, the off-diagonal block due to the influence of hydrodynamic on the left-hand 
side of the radiative equations can be stated 

1 f-u(E r +P r ) (E r +P r ) 



P 



A rh = -J - — J). (15) 



Since the jacobian matrix A is a block lower triangular matrix its spectrum can 
easily be computed and is composed of the eigenvalues of the hydrodynamic sub- 
system, i.e. u — a, u, u + a, where a is the speed of sound and of A_ + u and 
A + + u where A_ and A + are the eigenvalues of the Ml-radiation subsystem in the 
laboratory frame (cf. || and section 3.3). We remark that the nonconservative 
coupling terms (P r : Vu and F r .V u) appear only in the off-diagonal block A r h 
and thus do not influence the signal velocities of the whole system. 
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3. NUMERICAL TREATMENT 



The nonconservative form of the radiation-hydrodynamics system prevents us- 
ing classical Godunov type scheme without care. In this paper we focus on some 
applications such as radiative shocks or collapse of a gas cloud where various regime 
are encountered in the same simulations. 

Our numerical strategy is the following. The conservative system obtained 
by dropping the nonconservative coupling terms is dicretized by a Godunov-type 
scheme based on a Riemann solver and the nonconservative terms are discrctized 
using the same values as in this Riemann solver at the same time step. As we will 
see in the next sections the construction of the Riemann solver for the conserva- 
tive system is quite simplified by the lower triangular structure and merely reduces 
to the construction of two decoupled Riemann solvers for the hydrodynamic and 
radiation subsystems. 

The relaxation terms introduce an other kind of difficulty. For large value of 
the opacity they are stiff and an asymptotic analysis shows that the total energy 
is governed by a diffusion equation. This regime is commonly referred to as the 
equilibrium-diffusion limit ( pljj , JT3] ) . To capture this regime with the above system 
we must treat carefully the relaxation terms using an asymptotic preserving scheme 
following the ideas introduced in (||). This will be detailed in section 3.3.2. Finally 
equations (|l|), (^) and ( [Tof ) show that in a stationnary shock between two states 
where radiation is in equilibrium with material, the two states satisfy Rankinc- 
Hugoniot relations involving both the radiative and material quantities. We will 
show in sections 5.1 and 5.2 that our numerical scheme satisfies these relations even 
when there is large energy transfer between matter and radiation. 

In many astrophysical situations physical quantities vary over many orders of 
magnitude in a very steep way. To deal which such situations, the adaptive grid 
technic proposed by Dorfy and Drury @ (thereafter DD) is very well suited. We 
have therefore choosen to implement our radiative-hydrodynamical scheme on such 
a grid. However, the numerical schemes presented below can also be applied on a 
static eulerian grid which would be better suited for multi-dimcnsionnal simulations. 

Let us first start with some general remarks and notations. Space is discretized 
into computational zones centered on points r^. The grid interfaces are located at 
r i±i/2- The volume of a zone is given by: 

{|7r(r 3 ! — r 3 _ x ) in spherical geometry 
— r i-i = hi slab geometry (assuming a unit surface) 

Subscripts (i, i + i, i + 1...) always stand for the space coordinates, while super- 
scipts (n,n + i...) indicate time. Variables at time n + \ are defined by: 

X" + ^ = (1 - a)X? + aX? +1 a G [0, 1] 

for a — the scheme will be explicit and for a = 1 it will be fully implicit. We also 
introduce the following notations which are useful to have more compact equations: 

X [i±|] = (X i+| -^_x) and XW = (X n + l -X n ) 
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Subscript and superscript are always distributive: 

3.1. The grid 

In the following, we briefly summarize the adaptive grid method of DD. The 
grid concentration is defined by: 

L 

rii = — - - - 

where L is a typical length. If L is a constant, a uniform concentration corresponds 
to a uniform grid and if L — a uniform concentration corresponds to a logarithmic 
grid. This later case is useful for spherical geometry. 

The basic idea of DD is to demand that n is proportional to a resolution func- 
tion 1Z. This means that there will be a large number of gridpoints wherever the 
resolution function is large. The resolution function depends on the nature of the 
problems being solved, but it is most of the time a function of the magnitude of 
the gradient of the physical quantities in order to have a large number of gridpoints 
where there are steep features in the flow. 

In order for the grid to be computationally tractable it is also necessary to 
ensure that adjacent zones have roughly the same size and that the grid adaptation 
time is not too small compared to physical time scale. These two conditions are 
fulfilled using the following equation: 



7ii - a g (a g + l)(n t -i - 2th + n i+ i) (16) 
nl + ^W-hr 1 ) (17) 



where a g is the grid rigidity and T g the grid time-scale. Equation ( pL6|) ensures 
that the grid concentration does not vary from more than a factor (a g + l)/o; g 
between to adjacent zones and equation (|l7| ) set the time-scale on which the grid 
can readjust itself. 

The final form of the grid equation is then: 

& - 1 < 18 » 

This equation depends both on the grid coordinates and on the physical 
quantities entering the resolution function 1Z. 

3.2. The hydrodynamics 

In one dimension and taking into account the moving grid equations ([j])-(|3|) can 
be written: 

' d t pS + d r [pS{u-u g )] = 

dtpSu + d r [pSu(u — u g ) + PS] = Pd r S + S(£pF r + F) 
d t ES + d r [ES(u-u g ) + PSu] = -Snpc(a r T 4 - E r ) 

-S(fpF r + F).u 



(19) 
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where u g is the grid velocity and the surface S equals 1 in slab geometry and 4irr 2 
in spherical geometry. 

The left hand side of system ([l9]) is an hyperbolic system for the variables pS, 
pSu and SE. Now, by introducing the moment vector U = (p, pu, E) and the 
physical flux vector T{U) = (pu,p + pu 2 , u(E + P)), the hyperbolic part of system 
( [l9| ) can be written under the following form: 

d t (SU) + d r [S{T{U] - u g SU)} = 0. (20) 

As the physical flux vector is a degree one homogeneous vector valued function 
(ST(U) — T(SU)), if we set U = SU, the hyperbolic part (E0h can be rewritten as: 



d t U + d r [F{U) -UgUj = 0. (21) 
The system (^l| ) can be solved using an exact Riemann solver (Godunov scheme) : 

sK)r + (^ ) -^)^] =o (22) 

The state U* i is the solution on the line r/t — of the Riemann problem 
between the cell i and i + 1 involving the following initial conditions: 

(U(0,r) = S i+ ,U- H forr<0, 
\U(0,r) = S l+ iU+ H forr>0, 

where Sj , i is the surface associated with interface i + h. U~ i and U + , , arc 

the trace of U from the cell i and i + 1. For a first order scheme (no Van-Leer 
extrapolation) U7 i = £7, and C/^_i = f/i+i- (for details about the solution of the 
Riemann problem sec JL4|) 

We can note that we have used the solution of Riemann problem based on 
T(U) — u g U instead of F(U) flux. In fact, by the Galilean invariance of Eulcr 
equations, the solution of Riemann problem for T(Li) — u g U flux on r/t = is the 
same as the solution of Riemann problem for T(U) flux on r/t = u g . 

Unfortunately, the uniform flow is not a numerical solution of the scheme ( p2| ) . 
To satisfy this property, the mesh geometry must be solution of a balance equation 
called geometric conservation law (GCL): 

sH, -KU =0 ' (24) 

Of course, this relation is false for the spherical case. Moreover, if the product 
SiAri has been replaced by the volume of a cell Ai^, a new numerical scheme is 
obtained which satisfies a GCL balance law: 

1 / \ [n] / , 

aiH -KU =0 ' (25) 



I^ + ^-u.C |=0. (26) 
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The uniform flows are numerical solutions of the scheme (|26|). Hence, the whole 
numerical scheme about system ([l9]) can be written now under the following form: 

if^ + ^-^jpO (27) 



(d r S PAr + {^pF r + F)Av 



(28) 



^(^)' nI + (E*K-u 9 )+pv);+j ] 



(29) 



K/9c(a r T 4 - £7 r ) + (-pF r + F)u) Av) 

C ' J i 

3.3. Solving for the radiative transfer 

To simplify the presentation, we will first study the radiative transfer in a frozen 
medium and in a slab geometry. In this case, system (Q)-(@) gives 

d t E r + d x F r = ac(aT^ - E r ) (30) 
d t F r + c 2 d x P r = -acFr (31) 

where we have defined a = up. Using the Mi closure, the radiative pressure is 
given by : 

3 4 f 2 F 

P r = X E r = K where f = (32) 

5 + 2^4 -3/ 2 c:E r v ' 

The previous sytem is a hyperbolic one which can be writen under the following 
form: 

d t U + d x T(U) = S{U) (33) 

with 

»=(£:) **<>-(<&,) s ^'{" c{aT -i~F Er) ) («) 

IA is the moment vector, T the physical flux vector and S the source term vector. 
The Jacobian matrix J of this system is a function of / and c and is given by: 

The eigenvalues, A - , A + (A~ < A + ) of the matrix J can now be calculated: 



x'U)- Vx'(f) 2 - W(/) + 4x(/T ,„ v 

A = * c, (35a) 



A+ = x'(f) + Vx'(f) 2 - W(/) + 4 X (/T c (35b) 
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The eigenvalues, A~ and A + , of J are plotted in figure ffl- A]. At radiative equi- 
librium (when / = 0), the eigenvalues are equal to ±c/vo7 which is a direct con- 
sequence of the isotropy of underlying distribution function of the photons. On 
the other hand, in the case of extreme non-equilibrium, when / = 1 (resp. / = — 1 
), the eigenvalues are degenerated and are equals to c (resp. — c). This case cor- 
responds to the free-streaming limit where all the photons move towards the same 
direction. It is important to note that in the case of a constant Eddington factor 
(x = 1/3), the eigenvalues are also constant and equal to ±c/y/3, which means that 
a light front can not propagate at the speed of light, as it should. 




-0.5 0.0 0.5 -O.o 0.0 0.5 -0.5 O.O 0.5 

f f f 



FIG. 1 eigenvalues ofM\, — (full line) and — (dashed line) according to f — -jjj-> 
the normalized flux for e — 1 (A),e = 0.5 (B) and e = (C), 



3.3.1. Numerical approximation of P\ and M\ by the HLLE approximate Riemann 
solver 

The system (|30|)-(|3"l|) can be discretized using a Godunov method based on the 
HLLE approximate Rieamann solver (see Q). Hence, we need to find an estima- 
tion of the largest and smallest physical signal- velocities in the exact solution to 
the Rieamann problem. This estimation involves the eigenvalues of the Roe's lin- 
earization for the Mi model. The Mi Roe's matrix J ; , i between i and i + 1 states 
can be written as: 



3 i+i \c 2 ( Xa -X' m fa) C X ' r , 
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where 

• x(/i+i)-x(/i) 

/J t i Jt-\-± A,\J 't 1 ' A-VJlfJ-/ ' 



/< + /*+! .. _x(fi)+x(fi+i) _ I xt Xi-/- if /*^+i> 



[x'(/i) else. 

The eigenvalues, A~ x , A^ j (AT, x < A ^7 x ) of the matrix J; , i can now be cal- 
culated: 



\/xm 2 -4/ a x; n +4Xa 



c, (36a) 

- ' 2 Z 

a+ = y — a (36b) 



Once knowing the numerical signal velocities A* x , the numerical flux function 



± 

between cells i and i + 1 can be written under the following form: 

a+Fi-a-F i+1 + a+a-(U i+1 -Ui) , i a+ = max(X+ +1 , X+ ,,0), 

, i = : where < , 

+ 2 o+-a- a - = min(X i , A , , 0). 

The system (|3(])-(|3l]) can then be discretized in conservative form: 

Af + ^ = ^(a^ - (37) 

pn+l _ pn P"+5 - p n+ \ 

At + c 2 r;l+ \ x ^ = -v^cKP (38) 

3.3.2. Scheme modification to enforce asymptotic limit 

As we mentionned before, one of the important difficulties of radiative transfer 
is to be able to deal with very different regimes. It is therefore important to ensure 
that our scheme has the correct limiting behavior both in the free-streaming and 
diffusive regimes. 

In order to illustrate these points in a simple maner we will consider wave speeds 
equal to the speed of light (i.e. a + = — a~ = c) and a constant Eddington factor 
X = 1/3 (resp. 1) in the diffusive (resp. free-streaming) regime. The discretized 
equations then take the following form: 

^r.i ^r,i r r,i+l r r,i-l C , wre+i oc"+5 i c»+3 \ ™+5 I rri T? n +k\ 

At + 2K~x {E ^- 2Ei +E r,i-\) = °i 2 c(a r T m -E r ^) 

pn+l _ pn 771"+ 2 p"+2 

At X 2Ax 2Ax [ r ' i+1 r - i; CT * c ^ 

the modified equations actually solved by the numerical scheme arc therefore: 

d t E r + 9 x F r - °-AxdlE r = oc{aT^ n - E r ) (39) 
9 t F r + cV^r - ^d 2 x F r = -o-cF r (40) 
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Free-streaming regime 

Let us first consider the case of light freely propagating in vaccum (i.e. a — 0.) with 
all the photons going in the same direction (i.e. F = cE). The resulting equation 
for the energy of the radiation is then: 



d t E + cd x E r - —d 2 x E r = 0. (41) 

This equation is the exact advection equation with an added numerical diffusion 
term. The single mode solutions of the exact equation are of the form: E r — 
E r fi cos(27r/A(c£ — x)). For this mode and under CFL condition the variation of 
energy due to the diffusion term can be estimated: 



AE r _ Air 2 cAx AE r _ 4tt 2 cAt Ax / Ax \ 2 
~AT = T r '°~ ~E^ = ~2 A 2 VY) 

If the wave-length of the perturbation is properly sampled by the grid the diffusion 
term will always remain negligible. As could be expected, there will be some dif- 
fusion for very sharp energy profiles. In actual simulations, this diffusion will be 
further reduced by taking into account the proper waves speed and we will show in 
the next section that it remains acceptable. 

Diffusive regime 

Let us now examine the diffusive regime. In this case, the temporal variation of 
the radiative flux can be neglected. Equation (EOJ) then gives at first order: 



F r ~ -^-d x E r (42) 

OCT 



inserting the previous expression into equation (|39|) one gets 



dt E r -^-(l + ^)dlE r = ac(aT^-E r ) (43) 

The equation obtained for E r is indeed a diffusion equation. Furthermore, if 
the grid is sufficiently fine so that a Ax << 1, the diffusion coefficient is the right 
one. One the other hand, if the grid does not sample the photon mean free path 
{a Ax » 1), the diffusion coefficient is then dominated by a purely numerical 
term and it is therefore much larger than it should. With this simple analysis it is 
possible to explain some of our numerical results. For a coarse grid, the effective 
diffusion coefficient is too large and the thermal precursor is therefore longer. This 
wrong behavior in the diffusive regime is not a caracteristic of the HLLE scheme 
we have used. In any Godounov type method the divergence is discretized as the 
sum of a centered term and numerical diffusion term. This last term is absolutly 
necessary to ensure the stability of the scheme. Generally this term is proportionnal 
to Ax n , (n — 1, 2). When the mean free path of photon is much smaller than Ax, 
the numerical diffusion term can become dominant, as we have just seen. 

In order to have a scheme which reproduces correctly the diffusion limit, it is 
therefore necessary to explicitly take into account the smallness of the photon mean 
free path. We present below such a scheme. 



11 



We first write the system of equation (|3C|)-(|31j) in a dimcnsionless way. For 
that purpose we consider a characteristic length I and we defined the following 
quantities: 

e -_L t — — ~ — - F - — 
al al 2 I r ce 



the system ( |30| - |3lj ) can then be written: 

di E r + d £ F r =(°?^3± (44) 

diF r + \d iX Er = -\F r (45) 

The difficulty raised by the diffusive limit is highlighted in equation (flq). When 
e is very small, the hyperbolic system p4|)-(|4"5|) becomes very stiff because of the 
second term of equation fl45| ) (the jacobian condition number — > oo when e — > 0). 
An interesting way to overcome this difficulty is to transform the previous system 



in a non-stiff hyperbolic system with stiff source term j8||. Equation (45) can be 
written : 



dfFr + d^xEr 



: (F+(l-e 2 )d £X Er) 



(46) 



The new system to be integrated is now composed of equations ( |44| ) and (f46|). 
In order to integrate this system the hyperbolic terms (to the left of the equal sign) 
will be treated by the Godounov type method presented above, whereas the other 
terms (to the right of the equal sign) will be treated as source terms. Which means 
that they will be discretized in a purely centered manner. When e = 1 the sytem 
(^J)-(^) is identical to system (|30|)-(|3ll) and gives the correct transport regime. 
But when e — > 0, the source term in equation (||^) becomes dominant and gives the 
correct radiative flux in the diffusion limit. Furthermore, because of the change of 
F into F r , the numerical diffusion term in equation (^0|) has been multiplied by e 
and will therefore always remain small. 

The hyperbolic system to be solved is now: 



d~ t E r + dxF r = 
d- t F r + d iX E r = 



(47) 



As previously done, we will integrate the above system using a approximate 
HLLE Riemann solver. The Jacobian of the system (Et]) is: 



1 

x - fx' e x' 



with / = F r /cE r , x = X(f) and x' 



MS) 

df 



The eigenvalues, A + and A (A + > A ), of the Jacobian are given by: 
. ± ± V(ex') 2 +4(X-/X')) 
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It is worth noticing that in order to keep the previous eigenvalues real one must 
have (ex') 2 + 4(% — fx') > 0. Therefore, for a given value of / there is a minimum 
possible value (e m j n ) for the e parameter, (fx' — x) 1S an increasing function of 
/ which vanishes for f = fo = 2y3/5. Thus, if / < fo there is not constraint on 
the value of e (i.e. e m ,„ = 0) but if f > fo there is a strictly positive minimum 
value for e. The function e m in(f) IS plotted on figure (||). These limitations on 
e eventhough they have a mathematical origin, are physically justified since they 
imply that when the radiative flux is high one can not tend toward the diffusion 
scheme. Futhermore, we can deduce from system ( |47| ) and from the fact that the 
Mi model has an intrinsic flux limitation that / = F r /cE r < 1. It is therefore 
important to always keep e > / in order to compute correctly the radiative flux. 

In this case, we need also to find as in the previous subsection an estimation of 
the largest and smallest physical signal velocities in the exact solution to the Rie- 
mann problem. This estimation involves the eigenvalues of the Roe's linearization 
for modified Mi model (^). The Roe's matrix J i+ i for this model can be written 
under the following form: 

1 

(Xa -Xmfa) ^X'rr 

where we have used the same notation as in the previous section. The eigenvalues 
A. , i of this matrix are: 

4(Xa faXm) ) 

X *+h = 2 ' 

As in the previous case, one can give an estimation of the largest and smallest 
physical signal velocities: 

a + = max(Xf +1 , 1 , 0), a~ = min(X[ , X7 +1 , 0). (48) 
Using the HLLE solver gives the following expressions for the flux: 

~ = + j£l_ (Er . +1 Er . } (49 ) 



_ a+P rtl -a P r . l+1 5+5 - ~ 
' 2 a + — a a + — a 

Rewritting the system in the original variables, we get the following discretisa- 
tion for the Mi model: 



^(K +1 - + ^{F r ^ H - F r>i _i)"+i = ac(a r T 4 - E r )^ (51) 
±_ {F n + i _ F n h + l-(P r . H - = -° cF £ k ( 5 2) 

where we have defined: 



P r , i+h = e 2 P Ttl+h + (1 - € 2 ){P„ +1 + P„)/2 and F i± i = ceF r>i±i 
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FIG. 2 Minimum value of the e parameter in order to keep real the eigenvalues A + 
and A - plotted as a function of the reduce flux / = F r /cE r . 

The method is easily extended to the radiation subsystem written in the co- 
moving frame but where the nonconservative terms are dropped and with a moving 
grid. The sytem then writes 

d t E r + d x ((u - u g )E r ) + d x F r = ac(a r T 7 4 n - E r ), (53) 
d t F r + d x ((u - u g )F r ) + c 2 d x P r = -acFr, (54) 

The non-conservative terms are dropped for the present analysis but they must 
naturally be taken into account as will be detailed in the next section. The previous 
systeme can be written in a dimensionless way 

d~ t E r + d x {^E r + F r ) = (ar '" 2 " i H (55) 
d;F r + \d iX E r + -d £ (uF r ) = -\p ri (56) 

where 

u — u a 

u - 



c 

As for the system (0)-(§^) we transform the previous system in a non-stiff 
hyperbolic system with stiff source term 

d~ t E r + di{uE r + F r ) = 1 - E r ) - e(l - e)d i {uE r )) , (57) 

d~ t F r + di[ X Er + uF r ) = - 1 (P + (1 - e 2 )d xX E r + e(l - e)d £ (uF r )) , (58) 

whose eigenvalues are 

_ ± _ ± V(tx') 2 + 4(x - fx')) 
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The same analysis as above leads to the following scheme 

- E?)i + £(F* l+i - F^r+i = oc(a r T A - E r )^ (59) 



1 c 2 



At 

where 



_ + _ (P ; _ P * . )«+* = - ac F n r p (60) 



a + (u i+1 - u g )Ei - a (u i+ i - u g )E i+1 
F* i+l = F r . +| +e ±1 ' ,+ » *± (61) 



(1 



2 

.2; 



(u i+ i - u g )(E r , l+ i +E r ,i)), 



p :,+ l = e*P rAH + ^p-{P rii+1 + P rii ) (62) 
e a+((u l+ i - u g ))F„ - a~({u i+ i - Ug))F r>i+1 



c 2 5+ — a 



2c 2 



and 



max(A i t|_ 1 , AT^i , 0), a = min(A i , A i+i , 0). (63) 



In the previous discretization of the radiation subsystem, one need the value 
of e at each interface. In practical computations, where the opacity gradient can 
be large, we compute a value of e in each cell e = 1/aAx. The value of e at the 
interface is then taken to be the maximum value in the two adjacent cells. As we 
will see this method has given very satisfactory results even in the case of radiative 
shock with very strong opacity gradient. 

Let us also remark that if we take = X ± l , , then the HLLE solver coincides 
with the Roe solver for the radiation subsystem. Thus, if we use for the hydrody- 
namic subsystem a Roe solver, then the above scheme is exactly the Roe scheme 
for the radiation-hydrodynamics system where the non conservative coupling terms 
are dropped. Nevertheless using HLLE solver for the radiation insures positivity of 
the scheme. 



4. ACTUAL IMPLEMENTATION OF THE METHOD 



We present in this section how the full radiation-hydrodynamics system is dis- 
cretized and integrated. As we mentionned before, when using system (Q-^|) the 
treatment of hydrodynamics is independa nt of the radiation (apart form the source 
terms) and has been described in section (|3 . 2|) . For the radiative transfer in spher- 
ical coordinates the equations to be solved are the following: 



f d t (E r S) 
d t (F r S) 



d r (F r S) + d r ((u - u g )E r S) + P r SV.u 
npcS(a r T 4 - Er) + ^S(3P r - E r ) 



c 2 d r {P r S)+d r ((u- 
-npcF r S + c 2 P r d r S 



g )F r S) + F r Sd r u 

(3P r -E r ) n 2c 



(64) 
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The left-hand sides of the previous system contains the Mi-radiation subsys- 
tem and coupling terms due to the comoving frame. We treat the conservative 
terms using the Godunov-type method presented in the previous section. The non- 
conservative coupling terms, which, as we have seen previously, do not modify the 
eigenvalues of the physical system, are discretized in the following way: the velocity 
divergence is discretized using the value of the velocity at the interface given by 
the hydrodynamic solver. The two terms of the right-hand-side are local and are 
computed by integration over the volume of the cell. All the terms of the second 
equation are discretized in the same maner. And finally, all terms (given by the 
hydrodynamic or radiation solver, comoving and source terms) are computed at 
the same time so that all equations are fully coupled in an implicit scheme. This is 
summerized in the following tables. 



npcS \a r T 4 — E r ) 


P r SS7.u 


^S(3P r - E r ) 


(K P c{a r T 4 - E r )Av)" +i 


{Prt +i (OH] 


{f(3P r - E r )Av)^ 



npcF r S 


c Pf df S 


r 


F r Sd r u 


(KpcFrAv)^^ 


(c 2 P r S)^ 


( (3P ^c?Av) 


« 





In order to gain in accuracy, we actually use a classical second order MUSCL 
method (fill) to compute the flux. For the hydrodynamics minmod slopes are 
computed for p, u and P. For the radiative transfer, slopes are computed for E r 
and / in order to always satisfy the constraint |/| < 1. These slopes are then 
used to extrapolate the value of the different fields near the cell interface. These 
extrapolated values are then used as inputs for the Riemann solvers and for the 
upwind flux. 

Finally, as we will show in section 5.4, the CFL condition given by an explicit 
treatment of hydrodynamics is too much small in some applications and a full 



implicit treatment is needed. Then the non- linear system of equations (J1S|), (27)- 
( p9j ) and (Q) is solved iterativelly using a Raphson-Newton method. The jacobian 
matrix is computed using numerical derivatives. We use for the implicitation pa- 
rameter a — 0.55. 



5. TEST CASES 

We present in this section some test cases for the radiation-hydrodynamic scheme 
presented above. We start with some simple shock-tube problems and end with the 
collapse of a gas cloud leading to a stellar-like structure. All the shock-tube prob- 
lems are performed using a constant and uniform resolution (i.e. eulerian grid) in 
order to place emphasis on the properties of the radiation-hydrodynamics solver. 

5.1. Radiative front 

We first look at the propagation of a radiative front in vacuum. In the initial 
conditions the half-space x < is filled with radiation with T r = 10 4 and / = 1. For 
the half-space x > we have T r = 10 3 and / = 1/3. The front should propagate at 
the speed of light and remain sharp. The results are presented on figure (||) . We 



1G 



first see that, thanks to the Mi model, the front propagates at the proper speed. 
After a propagation of 1cm (ie 100 cells) we see that the front has suffered from 
significant diffusion and is spread among about 20 cells. But, when it propagates 
further the spatial extension of the front remains roughly stable. This stability of 
the front is due to the fact that for / = 1 the wave speed is equal to c while for 
/ = 1/3 it is equal to c/3. Therefore, contrary to a contact discontinuity, there is 
a slight compression of this radiative front. 

10000 
8000 

£ oooo 
E 

£ 4000 
2000 





-0.20 -0.10 0.00 0.10 0.20 

x/cm 

FIG. 3 The radiative temperature in the initial conditions (solid line), after a 
propagation of 1cm (ie 100 cells) (black cross), 2cm (dashed line) and 3cm (dotted 
line). Each curve is shifted along the x-axis of —cAt. Which means that since the 
front propagate at the speed of light the discontinuity should stay at x = 0. 

We now test our scheme with radiative shocks playing particular attention to 
the stability of the solution for different sampling of the photon mean free path. 

5.2. Radiative shock dominated by hydrodynamical pressure 

We examine in this section a stationary strong shock in a gas with unit molecular 
weight and a constant opacity of k = 4.0 10 4 err? j g. The upstream quantities for 
the gas are given by: 

p = 1.2 10~ 3 g/cm 3 T = 10 °K u= 200 km/s T r = T F r = 0. (65) 

As the upstream material is extremely opaque, any radiation incoming from the 
radiative front will be completely reabsorbed at some distance. Far from this area, 
the flow becomes homogeneous, the radiative flux vanishes (F r = 0), the radiation 
and the fluid temperatures are in equilibrium (T r — T). 

The downstream quantities are chosen so that the Rankine-Hugoniot jump con- 
ditions are satisfied for the mass ([j]), total momentum (||) and total energy (|l(J). 
Let us notice that in that case the Rankine-Hugoniot jump conditions for the com- 
plete coupled fluid radiative system (|l| - 1|) are not satisfied. As we introduce only 
three constraints on the downstream state, we need two extra relations to know 
this state (see Mihalas jll]]). The downstream material is also extremely opaque, 
hence, we still have F r — and T r — T . In fact, the gas given by the downstream 
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or upstream states act as a real non perfect gas where the material pressure P is 
replaced by the total pressure P + P r and where the material energy E is replaced 
by the total energy E + E r . In this test, the equation of state is totally dominated 
by the gaz (7 = 5/3) and radiation is only significant for energy transport. 

The photon mean free path in the upstream material is given by A = l/np ~ 
0.02cto. The results of the simulation for a grid with Aa; = 0.5A are given in figure 
(El). The hydrodynamical shock is stable and is sampled over 2 cells and there is 
a large region in front of the shock where the gas is preheated by radiation. In 
this thermal precursor the gas is heated and slowed down by the energy released 
in the shock. We have looked at the behavior of this precursor under changes 
in resolution. We can see (figure (|-A)) that even when the mean free path of 
photon is not sampled, the thermal precusor keeps the proper spatial extension. If 
the numerical scheme is not modified to take into account the diffusion limit, the 
extension of the precursor is much too large, as was expected from the analysis 
made in ( 3.3.2 ) . If the spatial sampling of the simulation is large compared to the 
length of the precursor, the temperature jump is captured over 2 or 3 cells which 
is roughly as accurate as the treatment of a shock. 



5.3. Shock dominated by radiative pressure 

We now examine the case of a shock where the dynamical influence of radiation 
is not limited to energy transport. In this test case, the energy in the downstream 
region is dominated by the radiative energy. The upstream conditions are: 



p = 1CT 3 g/cm 3 T = 10 5 °K u= 1000 km/s T r — T F r = 0. 

and the downstream conditions are given by the Rankinc-Hugoniot relations as 
detailed in the above section. As in the previous test we assume for the radiation 
that downstream states are at equilibrium, (i.e. E r = a r T^ and F r =0). 

The gas molecular weight is, as before, set to unity and the opacity is k = 
1 cm 2 J 1 g. The photon mean free path in the upstream region is therefore A = 
1/np — 10 3 cm. The results of this test case for Aa; = 0.5A are presented on figure 
. As in the previous case, the shock is captured in two cells and there is a large 
preheated region in front of the shock. In this region, energy is transfered from 
the gas to the radiation and the radiation energy density quickly dominates as the 
temperature rises. In the upstream region the equation of states is dominated by the 
gaz (7 = 5/3) and in the downstream region it is dominated by radiation (7 = 4/3). 
In this case the effective 7 in the downstream region is 4.2/3. Eventhough we do not 
explicitly treat the total (gaz and radiation) conservation law this energy transfer 
is treated correctly. As in the previous case, the shock is absolutly stable and the 
Rankine-Hugoniot relations are strictly verified. This example illustrates the fact 
that eventhough our hyperbolic analysis decouples matter and radiation, our scheme 
is valid even when there are large momentum and energy transfer between gas and 
radiation. In figure (0) we have plotted the temperature profile for different grid 
resolutions to illustrate the proper behavior of our scheme in both the streaming 
and diffusion regime. 



5.4. Collapse of a gas cloud 

Star formation is a very challenging problem, where radiation transport plays a 
dominant role. The modeling of the formation of a star introduces a wide range of 
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FIG. 4 Profiles of density (A), velocity (B), gas (solid line) and radiation (dashed 
line) temperatures (C) and radiative flux (D) for a strong shock with negligible 
radiative energy density. 



length and time scales. The mean density contrast between a collapsing cloud and 
the star reaching main sequence is about 10 20 . The initial cloud radius is about 
10 18 cm, while the typical radius of a solar-like star is less than 10 n cni. The com- 
pactness of the star embryo requires a spatial resolution of at least Ar/r ~ 1CP 4 
in the envelope outer layers, where the pressure gradients are the highests. This 
means that the Courant number corresponding to a dynamical time scale td = r j 'c Sl 
where c s is the sound speed, is about C = 10 4 . If the evolution of the system is 
to be followed over the Kelvin-Hclmotz time scale, tx-H = GM 2 /RL ~ 30Myr 
for a solar-like star, the Courant numbers we have to deal with is of the order of 
10 14 . These rough estimates demonstrate that the spatial resolution required for 
our problem demands an implicit numerical scheme. The process of star formation 
is dominated by the properties of radiative transfer in highly unsteady flows. The 
reason is that self-gravity alone can not lead to the formation of a compact object 
unless the energy released by gravitational contraction beeing radiated away, from 
the surface of the protostellar core. Energy transport in the core's interior, towards 
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FIG. 5 Profiles of gas temperature in the upstream flow preheated by the radiation. 
In plot A the grid spacing was Ax — 0.5A (solid line), 2.5A (dotted line), 5.5A 
(dashed line), 11A (solid line with squares). The solid line with triangle is the results 
of a simuation with Ax = 11A but where e was kept to unity. For the simulation 
of plot B the grid spacing was 550A, the thermal precursor is now sampled on 2-3 
grid points which is almost as accurate as the treatment of a shock. 

the photosphere, is mainly due to radiative transfer in an optically thick plasma, 
and convection. The photosphere is a narrow layer where the radiation transport 
regime abruptly turns from the diffusion regime to the free streaming regime, by 
which energy flowing from the core by diffusion and convection is radiated away. 
Furthermore, during the accretion phase, the photospere is surrounded by super- 
critical radiating accretion shock, where the kinetic energy of the infalling gas is 
essentially converted into outgoing radiation. Since most of the total luminosity 
of the protostar is produced within this sharp, optically thin accretion shock, the 
calculation of the structure of the radiative accretion shock is a crucial issue in the 
modeling of star formation. A correct evaluation of the entropy extracted by escap- 
ing radiation from the accreting material, deposited on the protostar surface, is also 
the key factor which determines the structure of the growing stellar core. First, the 
entropy profile determines the size of the quasistatic core, and hence determines 
the location of the protostar in the Hcrtzprung-Russell diagram by the end of the 
accretion phase; second the entropy gradient triggers the highly efficient convective 
energy transport, in regions where Vs < 0, where s is the specific entropy. 

We present in the following section an example of the calculation of the forma- 
tion of a primordial star, in a slightly bound region of an overall expanding universe. 
The simple model presented here includes only hydrodynamics, gravitation and, of 
course, radiative transfer, ft is nevertheless a very challenging test for radiative 
transfer. As we will see, the opacity gradients are very large and there is a steep 
transition from the diffusion to the free-streaming regime. 

We start with a 10 6 solar mass gas cloud. The velocities are those of the 
homogeneous expansion of the universe (i.e. u{r) = H r, where H is the Hubble 
constant). The average density is the critical density and the density profile is 
smooth with a small central overdensity. We start our simulation at a redshift of 
100 and the amplitude of the initial perturbation is adjusted so that the central 
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object forms at a redshift of 30. Initially the radiative energy density is that of 
the cosmic microwave background and the temperature of the gas is equal to the 
temperature of the radiation. 

In a first phase the entire cloud is in expansion. Then, at a redshift of 50 the 
central part of the cloud starts to collapse. The central density increases, but the 
temperature remains almost constant because the energy is evacuated by radiation 
since the photon mean free path is much larger than the typical scale of the forming 
structure. Next, at a redshift of about 30, the central region is dense enough to 
become opaque to radiation. The temperature and the pressure start to grow and 
become large enough to prevent the collapse in the central region. Finally the 
infalling matter creates an accretion shock on the central core. The structure of 
the flow when the central density has reached £.10 _8 <?/cm 3 is presented in figures 
©to®. 

The density profile is plotted is figure (g). The center is almost at the hydrostatic 
equilibrium and is surrounded by a slowly contracting envelope. The envelope ends 
at the accretion shock where the density jump is very important (about one order of 
magnitude) due to the fact that the gas radiates most of its energy, as is illusrated 
by figures (|h]) and (O). One can notice that the density contrast between the 
center and the outer accretion flow is of more than 18 decades and that the length 
scale vary of about 10 decades. 

The velocity profile (figures ([)])) shows that the outer region of the cloud is 
still expanding, following the global expension of the universe. On the contrary, 
the inner region is gravitationally bounded and is collapsing. The accretion shock 
is clearly visible and is sampled over only one computational grid point. Due to 
the moving grid, the central part and the region of the shock have a high spatial 



resolution. The actual resolution profile can be seen on figure (13). In the central 
part of the cloud, the radiation is in the diffusion limit and the luminosity slowly 
increase (figure (JTc|)) . At the shock there is a very sharp increase in luminosity due 
to the fact that the gas is heated at the shock and radiates most of its energy. The 
luminosity is then constant except for the outer region where the luminosity of the 
cosmic microwave background starts to dominate. As is illustrated by figure (|ll|), 
the radiation is in the diffusion regime in the central part of the cloud where the 
photon mean free path is much smaller than the radius while it is clearly in the 
free-streaming regime in the outer part. The transition between these two regimes 
occurs in a very sharp way at the shock. Figure (|lj) shows a zoom of the gas and 
radiation temperatures profiles near the shock. The spike in the gas temperature 
is a characteristic feature of supercritical radiative shock and a very high spatial 
resolution is needed to resolve it. Finally, the profile of the Eddington factor is 
plotted on figure ([l2]). The Eddington factor is equal to one third in the central 
object, which confirm the fact that the diffusion regime is valid in this region. Then, 
there is a small jump at the shock and the Eddington factor grows as one goes away 
from the surface of the star. At large distance, the light coming from the star is 
blured into the cosmic microwave background and the Eddington factor decreases 
back to one third. These large variations of both the Eddington factor and the 
opacity illustrate the interest of the M\ model for such radiative shock. 

6. CONCLUSION 

We have presented in this paper a new numerical model for radiative transfer 
which gives satisfactory results in a wide range of physical conditions including both 
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the diffusion and the free-streaming regimes where the anisotropy of the photons 
distribution function can be large. This method is therefore well suited to mod- 
eled some laser experiments and many astrophysical flows. Works are in progress 
to incorporate in the model a more detailed description of the physical processes 
(multigroup approach, ionization, ....) and to further develop this approach in a 
multidimcnsionnal framework. 
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FIG. 6 Profiles of density (A), gas (solid line) and radiation (dashed line) tempera- 
tures (B), radiative flux (in arbitrary units, solid line) and reduce flux (dashed line) 
(C) and gas (solid line) and radiative (dashed line) pressures (D) for a strong shock 
where the downstream energy is dominated by radiative pressure. The squares for 
the gas density and pressure show that eventhough the radiative pressure dominates 
the hydrodynamical shock is still properly treated. 
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FIG. 7 Profiles of temperature in the upstream flow preheated by the radiation. 
The grid spacing was Ax = 0.5A (solid line), 35A (dotted line), 200A (squares). 
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FIG. 8 Density profile of the forming star. 
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FIG. 9 Velocity profile of the forming star. 
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FIG. 10 luminosity profile of the forming star. 
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FIG. 11 This plot represents the ratio of the photon mean free path over the 
radius as a function of the radius. This ratio is very small in the center where the 
coupling between matter and radiation is very strong, and it is very large in the 
outer region where the gas is transparent. This parameter varies over 10 orders of 
magnitude, with a sharp discontinuity at the shock, where the transition between 
the diffusion and the free streaming regime occurs. 
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FIG. 12 Profil of the Eddington factor. 
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FIG. 13 Profil of the grid resolution. The grid resolution is very high at the 
shock, where all the gradients are very large. With the high resolution it is possible 
to see the temperature spike at the shock (see figure (fj~4|)). The resolution is also 
important in the center to properly sample the pressure gradient and near the peak 
of the Eddington factor. 
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FIG. 14 Temperature profiles at the shock. The gas is hcaten by the shock and 
then very quickly radiates its internal energy. This results in the large spike in the 
gas temperature at the shock which is associated to the jump in luminosity. On 
the contrary, the radiation temperature varies smoothy through the shock, it only 
changes its slope. 
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